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The behaviour of a system of ion channels formed across the cell mem- 
brane is presented. The infinite number of channels with an infinite cou- 
pling is introduced first as a reference point for the detailed derivation of 
the thermal-equilibrium probability distribution and the classification of 
the phase transitions. Fluctuations in a finite system are discussed next. 
We propose a new, step-like model of the ion channel switching, for which 
we provide the analytical results. The relation of this model to experiment 
is also provided. Finally, the master equation for the finite channel-number 
membrane is analysed numerically with the help of an exact- diagonalization 
technique. In particular, the decay-time of a metastable solution is esti- 
mated. The results do not agree with those obtained perturbationally, the 
difference is explained by the proposed frozen- diffusion approach. 

PACS numbers: 87.10.-He, 87.16.Uv, 05.40.-a 



1. Introduction 

All living cells are effectively separated from the environment by the 
lipids that constitute the cell membrane. The lipid bilayer formation and 
its dynamics have been recently objects of an intensive molecular dynamics 
studies [1]. However, ions such as sodium Na"*" or potassium K"*" can interact 
with proteins in the cell membrane that transport the ions across the mem- 
brane. Ion channels are one of the three types of proteins involved in the 
ion transport; the others are the pumps and the carriers. Approximately, 
an equal portion of the ionic current passes through each of these routes [2]. 

Classic experiments using the patch clamp technique showed an inter- 
esting stochastic behaviour of the ion channels that is controlled by the 
voltage across the membrane [3]. For example, this behaviour is crucial for 
the phenomena of neural excitability [4]. The fractal model of ion channel 
kinetics was formulated [5] and widely discussed in relation to the Markov 
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model [6]. Recently, Kramers' diffusion theory has been considered as an 
alternative to those models [7]. One should also mention a significant im- 
provement in the data analysis [8] namely, a new method was developed 
for the identification of deterministic structure in the time series, such as 
obtained from the patch clamp experiments, and successfully applied to the 
hearth disease and the epilepsy research. 

In all above experimental and theoretical analyses, the kinetics of the 
single ion-channel was investigated. In this paper we discuss collective ef- 
fects in both infinite-size system (Sec. 3), and realistic system (Sec. 4-6) of a 
few-hundred coupled ion channels along the cell membrane. We discuss the 
so-called step-like empirical model of the ion-channel switching afi^ected by 
the voltage across the membrane that is based on that proposed by Hodgkin 
and Huxley [4], but is much simpler and approximately rationalises the ex- 
perimental data for the potassium and the sodium channels [9] . 



2. Model and the physical system 

In this section we discus a simple model of ion channel switching and 
apply it to the system of m-channels across the membrane (see Fig. 1). 

The ion channels are formed by specific membrane proteins which un- 
dergo a spontaneous, voltage-sensing conformational transitions between the 
open and the closed states [10]. Hereinafter, we will concentrate on isolated 
assembly of m ion channels of one type, i.e. potassium K^. Let us suppose 
the simplest 2-state Markov model of the single channel switching (this is 
not a strict supposition, because such a model is equivalent to the m+1 - 
state Markov model of the whole membrane, cf. Sec. 4). In this model such 
a channel is characterised by two microscopic parameters determined by 
the voltage Vm across the membrane: the inverse average open- and closed 
times {Xopen = Aopen(Ki), and Xdosed = KlosediVm), respectively). Physi- 
cally, parameters Xopen and Xdosed are measured indirectly in a time scries 
obtained from patch clamp experiment for single channel [3], whereas Vm 
is the driving parameter controlled in such an experiment. Subsequently, 
the expectation value of the fraction of opened channels (uopen) obey the 
two-state, "open-closing" dynamics [4]: 

{iT'open} — Xopen (1 {l^open}) Xdosed (j^operi) • (1) 

For the sake of convenience, one can define the decay constant A = Xopen + 
Xdosed (that provides natural timescale of the processes) , and the probability 
of single channel open Popen = Xopen/ X, which allow us to rewrite eq. (1) in 
form: ^ ^ 

^ {iT'open) ~ Popen {f^open} ■ (2) 
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In further considerations, parameter A = A {Vm) is regarded as constant, 
and used to define the dimensionless time: 

T = At. (3) 

One can verify this approximation for a broad range of voltage across the 
membrane Vm by looking at the experimental data (see Table 1). With this 
simplification, only one microscopic parameter Popen {Vm) remains. 

Our starting assumptions, concerning the electrical properties of the 
membrane (see Fig. 1), are as follows: (i) the conductance of the membrane 
is directly proportional to the number of open channels, i.e. its resistivity 




Fig. 1. The patch clamp experiment version for a m-channel membrane and the 
corresponding electric circuit. The membrane parameters are: Ip - the constant 
current of the ionic pumps, Rq - the resistivity for the ionic diffusion across the 
membrane, C - the membrane capacity, and Rm - the m-channel system resistivity 
when all the channels open. For further considerations we suppose, that the system 
is connected to the external voltage V and resistivity R. Below: the effective circuit 
for long-timescale processes in the membrane {t RaC). The effective external 
voltage and resistivity R^s are defined in the main text. 
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Table 1. Microscopic parameters of single potassium channel from human ocular 
epithelial cells (after Twohig, Ref. [9]). 



Membrane 


Ave. open 


Ave. closed 


Time 


Prob. of 


potential 


time 


time 


unit 


open chan. 


Vm [mV] 


1/Aopen [ms] 


1/Adosed [ms] 


1/A [ms] 


Popen 


-20 


14.57 


1.58 


1.425 


0.098 


-30 


11.84 


2.04 


1.740 


0.147 


-40 


5.18 


3.62 


2.131 


0.411 


-50 


5.42 


5.42 


2.710 


0.500 


-60 


3.94 


8.06 


2.646 


0.672 


-70 


2.72 


15.46 


2.313 


0.850 


-80 


1.77 


39.37 


1.694 


0.957 


-90 


1.71 


60.32 


1.663 


0.973 



is equal to Rm/nopen, where Rm is the resistivity when all the channels 
arc open, (ii) the total current of the ionic pumps Ip is constant, {in) 
the voltage-induced ionic diffusion across the membrane is described by 
the linear resistivity Rq, (iv) the relevant time-scale is that of long-time 
processes t » RqC, where C is the membrane capacity, and {v) whole the 
system {i)-{iv) is coupled to the constant voltage V source with internal 
resistivity R. 

After introducing all above suppositions, a direct application of the Kir- 
choff 's laws to the electric circuit shown in Fig. 1 proves its equivalence to 
the effective circuit, which contains only an assembly of switching channels 
coupled to the external voltage V^g and the resistivity Res: 

V-RIp „ R 

= 1-R/Rg' = 1-R/Rg- 

The calculation of the relation between the voltage across the membrane 
Vm and the fraction of open channels Uopen for the effective electric circuit 
is straightforward and yields 

J- "T ii'open^^S / 

This allows us to consider the probability of a single channel open as a 
function Popen = Popen (nopen), sincc Vm is & function of Uopen (5). Moreover, 
instead of applying phenomenological models of single channel switching 
Popen (Kn), such that proposed by Hodgkin and Huxley [4], we can now 
consider indirectly a model function Popen {nopen) , describing the assembly 
of coupled ion channels themselves. 
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Namely, for a qualitative analysis of the membrane, one can select any 
model function Popen (nopen) that retains principal properties of physical 
system, such as negative feedback (which means that channel remains closed 
for high polarising voltage Vm and open for low, see [9]), and have convenient 
analytical properties. In Sec. 4-6 we apply the step-like function 



Popen iP'open) — ^1 + 6 



B(no 



) 



(6) 



where B and no are adjustable parameters, which allow us to perform some 
integrations analytically. A brief comparison of the model of channel switch- 
ing defined by the function (6) with some of available experimental data [9] 
are provided in Fig. 2. Strictly speaking, the relation inverse to (5) per- 
mitted us to write the probability of channel open as a function of voltage 
again, i.e. 

Popen(Kn) = (l + e^+''°/'''")"\ (7) 

where A = B{no + Rm/Res) and Vq = —BVesRm/Refi arc empirical param- 
eters, constant for a particular type of ion channel (i.e. potassium channel 
for human ocular epithelial cells, see Fig. 2). Empirical parameters A and 
Vq provides the relation between model parameters B, no and the physical 
driving parameters Ves, Rm/ReS (see Fig. 1): 



Bno-A' 



R^/Res = A/B - no. 



(8) 
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Fig. 2. Comparison between the step-like model of ion channel switching for the 
fitted empirical parameters A = — 3.85±0.78 and Vq = — 142±30 mV (see main text 
for a detailed description) and the experimental data for human ocular epithelial 
cells listed in Table 1. 
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In effect, the selection of the probability of channel open in the form (6) 
constitutes the step-like model, which, as we just shown, rationalises the 
experimental data of Ref. [9] . 

3. Deterministic system: long-range coupled infinite membrane 

In this section we consider the limiting case of the membrane with infi- 
nite number of channels {m — > oo), and far from a critical point [11]. For 
such a system one can suppose that the probability distribution is narrow 
and, subsequently, the fraction of open channels Uopen is close to its aver- 
age value Uopen ~ {iT-open)- With this assumption (which will be verified in 
the next section) , we can write down an approximate version of Eq. (2) as 
follows ^ 

"Ij^'^open ~ Popen ij^open) "^openi (9) 

where r is a dimensionless time (3). Equation (9) becomes exact in limit 
m ^ oo (except in the case when the system undergoes the 2-nd order 
phase transition), as a mean- field solution of the systems with the infinite 
coupling [11]. 

For all the models of ion channel switching [9], qualitatively similar to 
that defined by (6) , eq. (9) could have from one to three stationary solutions 
iT-openiT) = thsit ohej the relation (cf. Fig. 3) 

Popen (n*) = ra*. (10) 

For each of these solutions, one can consider its stability against small fluc- 
tuations Uopen = -|- 6u, by making a proper linearization of the relation 




Fraction of open channels n, 



Fig. 3. Stable and unstable stationary solutions n* = Popen {n*} of the dynamical 
equation for the infinite membrane (see main text). The step-like model of ion 
channel switching is applied. 
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(10): 

Popen (nopen) = Popen + Sn) = + p'ju + O (Sn^'j , (11) 

where p'^ is the first derivative of the probability as a function of the fraction 
of open channels: p'^ — Popen ('^*)- The equation (11), after substituting to 
the dynamical eq. (9) provides us with an approximate result (for Sn « 1): 

5n{r) = 5n{0)e-^^-P'*^^. (12) 

This result simply means, that the first derivative of the probability p'^ in 
the stationary point defined by the relation (10), allows us to identify 
stable (for < 1) and unstable (for p'^ < 1) stationary solutions of the 
dynamical eq. (9) (see Fig. 3) . The threshold value p'^, = 1 will be shown to 
correspond with the critical point in terms of 2-nd order phase transition in 
the next section. 

The approach presented above provides us two equivalent stable solu- 
tions (see Fig. 3). In the next section we will start the considerations for 
finite number of channels m, and will verify (in Sec. 5) one of these stable 
points to be actually metastable, with a finite decay time growing exponen- 
tially with m. 



4. Stationary probability distribution for a finite membrane 

Let us consider the finite, m ion channel cell membrane, and the evo- 
lution of its probability distribution P = {Pk}, (where k = 0, . . . ,m is the 
number of open channels, and the normalization condition J^k -Pfc = 1 is 
satisfied) in time. All the probability flows between fc-th point and two 
neighbouring points are shown symbolically on Fig. 4. 

In the limit of the stationary probability distribution P^'^^ = |P^''^| one 
can write down the detailed-balance equation between k — 1 and k points 
(cf. Fig. 4) in the form 

{m-k + l)pk-iPl% - k{l - Pk)Pt''^ = (13) 

for each k = 1, . . . ,m, where pk is the probability of open channel when k 
channels are open, pk = Popenik/m) (see Sec. 2). For non-zero probability 
distribution (this seems to be natural supposition for physical systems and 
will be verified a posteriori) one can rewrite Eq. (13) in the form 

Pf = ^ -'^+ 'P'^-^ p(% (14) 
" k l-Pk 

for k = l,...,m, that determines the stationary probability distribution 
p(o) = |Pj^'^'*| except of the multiplicative constant Pq'\ which has to be 
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adjusted to conserve normalization Y^j^ Pj^ = 1. Eq. (14) provides the pow- 
erful tool for estimating numerically the stationary probability distribution 
for m-channel system with given model of ion channel switching Popen{n). 
Moreover, in Appendix A, using Eq. (14), we discuss small fluctuations 
around the stationary solutions (10) for a general model Popenin-)- 

The numerical results, presented in Fig. 5, as well as an approximate 
formula for small fluctuations (A. 8), derived in App. A, suggest the approx- 
imate analytical way of solving detail balance equations for large m. In this 
limit, one can look for universal function f{k/m), such that 

Po(n) = e"*(^(")+^'"). (15) 

In Eq. (15) we define a quasi- continuous variables: fraction of open channels 
n = k/m, and the stationary probability distribution function Poin), defined 

by relation pjf^ = Po{k/m)An, where An = l/m. Cm is an adjustable 
parameter, so that function Po{n) (15) satisfy the normalization condition 

Jq dnPo{n) = 1. One can roughly estimate the value of Cm by making 
Gaussian approximation around the global maximum of /: f{n) « fmax — 
fmaxin - n^ax)V2, obtaining: 

1 In ?77 

Cm = -fmax + 7^ + 0{l/m). (16) 

Vk-iPk-iAr -VkPkAr 



-UkPkAr Uk+iPk+iAr 




Uk - k{l - Pk) 
Vk = {m- k)pk 
Pk = Popen{k/m) 



k-1 k k+1 



Fig. 4. Time evolution of the probability distribution for finite m-channel mem- 
brane toward the stationary state |-Pfe°''|- Index k = 0, . . . ,m is the number of 
open channels, the amounts of probability flowing between each pair of the points 
in the dimensionless time interval At are shown. The coefHcients Uk, Vk and pk in 
probability flows are defined above. 
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Then, making an expansion 



mf{n — An) = m 



f(n)- ^An + 0(An^) 
an V / 



= mf{n)-^ + 0{l/m), 
an 

(17) 



and substituting it into the detail-balance equations (14), we obtain: 

-4L _ n 1 - Popciin) 



1-n Popen{n) 



+ 0(l/m) 



(18) 



Eq. (18) in turn, provides the approximate differential equation for func- 
tion /(n) (exact in the hmit m — ^ oo): 



^ ^]^(^~^ Popenin) 

dn 



n 1 - Popenin) 



(19) 



-0.06 




I(n) + L.1000 



0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 
Frac. of open channels n 



Fig. 5. Exact numerical solutions of the detailed-balance equations for different 
m compared with an analytical solution of corresponding differential equation for 
m — > oo (see main text). The quasi-continuous probability distribution function 
is defined as P{n) = mPk, where n = k/m. Above: symmetric step-like model 
of ion channel switching {B = 5, no = 0.5), below: an asymmetric case {B = 5, 
no = 0.49) 
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In particular, when looking for extrema of the function /(n) (the condi- 
tion df/dn = 0), which are equivalent to the extrema of the stationary 
probability distribution function Po{n) (15), Eq. (19) reduces to the con- 
dition for stationary solutions of the deterministic equation of motion (9) 
JT-* = Popen{i^*) (see Sec. 3). This provides the link between these two ap- 
proaches. 

Eq. (19) can be solved analytically for the step-like model of ion channel 
switching Popenip^) defined by Fermi-Dirac function (6): 

f{n) = —n In n — (1 — n) ln(l-n) + -(n-no)2. (20) 

Plots of the function (20) for sample parameters B and no are shown in 
Fig. 6. According to the formal correspondence of Eq. (15) with the prob- 




B^ = 5, no= 0.50 -----y- 

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 
Frac. of open channels n 

Fig. 6. Free energy per channel in the step-like model of ion channel switching 
for different parameters B and rip. Top panel: 2-nd order phase transition for 
symmetric case (no = 0.5) and varying B, Bottom: the phase transition for fixed 
B = 5 and varying no {note the presence of the secondary minimum on the left, 
corresponding to a metastable state). 
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ability distribution for Gibbs canonical ensemble [11] 

P, = = e-^in,B,no)/e^ (21) 

where is the sum of states above all (jjfj configurations for given number 

of open channels k, O is the reduced temperature, and J^{n, B,no) is the 
free energy of the system with given fraction of open channels n, one can 
formally define the mean- field Ginzhurg- Landau [11] free energy functional 
for m — ^ GO limit (and subsequently {n))'- 

J^in,B,no) = -mfin). (22) 

Here we will concentrate only on the isothermal ensemble (0 = 1), [12]. The 
definition (21) allows us to identify the candidates for 1-st and 2-nd order 
phase transitions (cf. Fig.6) in terms of Ginzburg-Landau theory [11]. How- 
ever, when the parameter B is constant and uq varies, the configurational 
entropy of the system 5conf = In does not change when the fraction of 
open channels n switch between two equivalent minima at the point n = 0.5, 
so this process also has to be considered as the 2-nd order phase transitions 
within the Ehrenfest classification scheme. The remarkable feature of this 
case is possible existence of a metastable state (cf. Fig.6), that motivated 
us to study the dynamics of the system in the next section. 



5. Dynamics of the probability distribution via an exact 

diagonalization 

In this section wc derive the system of dynamical equations (the master 
equation), which describe time evolution of the probability distribution of 
the number of open channels for the m-channel membrane. We will also 
present results obtained when solving this system by using the numerical 
exact diagonalization technique. These results are substantial for further 
considerations, i.e. they are used to justify the approximate approach pro- 
posed in Sec. 6 for large systems. 

Let us consider a discrete probability distribution for the m-channel 
membrane (Pk), where A: = 0, . . . ,m is the number of open channels, and 
the normalization condition J2k -Pfc = 1 is satisfied (see Sec. 4) . Taking into 
account probability flows between the k-th point and its two neighbours 
k — 1, k + 1 (see Fig. 4) one can write down the system of linear, ordinary 
differential equations: 

T m 

-^Pk = -Y.AkiPi, (23) 
"''^ 1=0 
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where matrix [A^i] has a tridiagonal form A^i = {uk + Vk)Ski — Uk+i6k+i,i — 
Vk-i6k-i^i, 5ki is the Kronecker delta, Uk = k{l — Pk), Vk = {m — k)pk (see 
also Fig. 4), and the probabihtics = Popenik/m) are given by the model 
of ion channel switching (cf. Sec. 2). 

The system (23) can be simplified by the rescaling Lk = Pk/o-k, where 
coefficients (a^) are defined for = 1, . . . , m by the recursive formula: 

ak = .l^=(^^^^'-^Y\ (24) 
\ Uk \ k Pk J 

and aQ remains a free multiplicative parameter. However, comparing the 

formula (24) and the detailed balance Eq. (14), one can choose = yPjF^ 
loi k = 0, . . . ,m (cf. Sec. 4), and adjust ao appropriately. That form will 
be suitable for e.g. perturbation approach (see Sec. 6). In terms of variables 
(Lk), the linear system (23) takes the form: 

7 m 

Lk = -J2BkiLi, (25) 

"''^ 1=0 

where [Bki] is the symmetric, tridiagonal matrix with diagonal elements 
Bkk = Uk + Vk, and off-diagonal Bk-i^k = -akUk/ak~i = -yJVk^iUk- For 
particular model of ion channel switching, in our case the step-like model (6), 
the matrix [Bki] can be diagonalized with standard numerical packages [13], 
up to m « 1000. Presented technique allows us to examine time intervals 
as long as r Ri lO^'^ (3), unreachable by an indirect numerical methods of 
solving of the system (23) (e.g. with Runge-Kutta method). These features 
have substantial meaning for searching metastable states for a large system 
that will be discussed further in this section. 

The formal solution of the system (23) is given by the linear combination: 

m 

Pfe(r) = afe^q4')e-'^ (26) 
/=o 

where lji are the eigenvalues of [Bki] in ascending order ljq < ... < ujm, 
{L^'^} are the corresponding, orthonormal eigenvectors, L*-^-* = k, I = 

0, . . . ,m (so that J2k ^k^^k^^ ~ ^in), and coefficients q are defined by de- 
composition of the initial probability distribution {-Pfe(r = 0)} in the basis 
{LW}: 

m 

ci = ^4^Pk{0)/ak. (27) 

k=0 
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One can find the lowest eigenvalue of the matrix [Bki] to be equal ujq = 0, 

by substituting L^*^-* = \fp^ , for /c = 0, . . . , m, where P^^^ = is the 

stationary probability distribution obtained in Sec. 4 from detailed-balance 
equation (13). In other words, the solution (26) of the system (23) relaxes to 
the stationary distribution |-P^'^^|, for r — > oo. In all the numerical analysis 
below, we examine the details of this relaxation, particularly we search the 
possible metastable state associated with the secondary maximum of the 
stationary probability distribution {Pf)} (see Fig. 5). 

We have applied the step-like model of ion channel switching (6) in 
symmetric {B = 5, no = 0.5) as well as asymmetric {B = 5, uq = 0.49) 
case. The initial state (Pfc(O)) was chosen to be the 50% mixture of the 

stationary probability distribution (j^^^''^ and 50% of narrow-peaked distri- 
bution, concentrated in secondary maximum of the stationary distribution 
(see. Fig. 5). Two characteristics periods have been identified in the time 
evolution of m = 400 channels membrane, in both symmetric and asym- 
metric case: (i) a fast diffusion from the initial conditions to metastable 
solution (see Fig. 7), and (ii) a slow relaxation to the stationary solution. 
Both of these features arc shown in Fig. 7. Our quantitative observations 
for the membranes with different numbers of channels m can be briefly re- 
ported as follows: the time of diffusion (i) was found to be approximately 
constant for m = 100-^400, whereas the life time of the metastable solution 
(m) grows rapidly with m. Moreover, for large systems, the shape of the 
probability distribution remains exactly constant for a long period of time, 

1. e. T = 10^ 10^ for m = 400, that correspond to the real time approx. 
0.2^20 s (see Table 1). 

According to these results, one can expect a broad window in the spec- 
trum of eigenvalues (0;^) of the dynamical system (23). Spectra presented in 
Fig. 9 confirm clearly this hypothesis, for both the symmetric and the asym- 
metric cases. Moreover, first excited-state eigenvalue uji increases exponen- 
tially with the number of channels m, (except of the numerical round-off 
errors, which become significant for coi < 10~^^), whereas all the other 
excited-state eigenvalues lie in approximately constant (for different m) 
band bordered by the eigenvalues u;2 and tOm- This simple picture proves, 
that for the m-s of order few hundreds and larger, the dynamics of the prob- 
ability distribution (P^) is completely determined by first eigenvalue loi and 

associated eigenvector 

p(i) = (^pwy 

for times longer than r ^ 100 (that 
correspond to the real time approx. 0.2 s). 

A further numerical evidence for this statement is provided with Table 

2, as well as Fig. 10. and 11. Let us define the trapping times T^^) ^^'^ 
Tsh as a position of the last inflection points (before the system reaches 
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Table 2. Brief comparison of the decay-time defined by first excited eigenvalue 
l/u)\ with trapping times T(„) and Tsh (see Fig. 10 and 11 for definitions). Sample 
results for symmetric {B = 5, no = 0.5) and asymmetric {B = 5, no = 0.49) 
step-like models of ion channel switching are provided. 



no 


m 


- log UJi 




^ogTsh 


0.5 


100 


2.748 


2.540 


2.392 




200 


4.293 


4.106 


3.968 




400 


7.397 


7.230 


7.091 


0.49 


100 


2.364 


2.195 


2.110 




200 


3.236 


3.057 


2.948 




400 


4.942 


4.768 


4.653 




Fig. 7. Evolution of the initial probability distribution for m = 400 channels mem- 
brane (see main text for details) to the metastable solution, obtained from the 
numerical exact diagonalization technique. The stationary solution for r ^ oo is 
also presented. The quasi-continuous probability distribution function is defined 
as P{n) = mPk, where n = k/m. Top panel: symmetric step- like model of ion 
channel switching with no = 0.5 and B = 5. Bottom: an asymmetric case with 
no = 0.49 and B = 5. 
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an equilibrium) of average fraction of open channels (n) and the Shannon 
entropy Ssh, respectively (see Fig. 10 and 11 for details). By looking at 
Table 2, one can find the trapping times T^„-) and Tg^ to be in the same order 
as corresponding decay times defined by first excited eigenvalue l/wi, i.e. 
Tsh ~ 10~'''^/ti;i « l/2a;i, where the factor 1/2 could be simple explained by 
representing the probability distribution in the form: P{n) = Po{n) + S{n), 
and calculating the Shanon entropy in this manner: 

Ssh= [ dnP{n) In P{n) = So -I f dn-^ + 0{S^) So + const e'^^^^, 
J 2 J Po{n) 

(28) 

where approximation works for r S> l/a;2 and So is the value of Shanon 

entropy for the stationary distribution Pq (n) . One remarkable feature of the 
Shanon entropy Ssh for the membrane is that its value for the stationary 



T = 10^- 10^ . 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 
Frac. of open channels n 

Fig. 8. The metastable probability distribution for m = 400 channels membrane ob- 
tained from numerical exact diagonalization technique, (datapoints does not change 
visibly in given time intervals) and its relaxation to the stationary solution for 
r — » oo (also presented). The quasi-continuous probability distribution function 
is defined as P{n) = mPk, where n = k/m. Above: symmetric step-like model 
of ion channel switching {B = 5, no = 0.5), below: an asymmetric case {B = 5, 
no = 0.49). 
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distribution is not always maximal, for same cases entropy Ssh can decrease 
as a function of time (see e.g. Fig. 11). This is caused by the fact, that we 
are studying a dissipative system, connected to the environment (see Sec. 2), 
whereas only the total entropy of the system together with its environment 
grows in time, not necessarily the fractional entropy of the membrane (i.e. 
Shanon entropy). 

Observations discussed in this section, particularly the simple structure 
of the eigenvalue spectrum (presented on Fig. 9), gave us the motivation 
for some approximate approach for the first-excited eigenvalue uji, that is 
discussed in the next section. 
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Fig. 9. Spectrum of eigenvalues obtained from the numerical exact diagonalization 
approach, as a function of the number of channels m. The first (wi), the second (^2) 
and the last (wm) excited eigenvalues are shown, all the others are lying between 
u)2 and uim (except from loq = 0). To make interpretation in terms of decay time 
easier, — logwfe are presented. Above: symmetric step- like model of ion channel 
switching {B = 5, no = 0.5), below: an asymmetric case {B = 5, no = 0.49). 
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6. Perturbative analysis and discussion 

Let us consider the membrane with large number of channels (m = 
200^ 1000, e.g.) and a metastable solution for the fraction of open channels 
um clearly separated from the stable ns (in other words: secondary and 
primary maximum of the stationary probability distribution, see i.e. Fig. 5.) 
by an unstable point nu (and corresponding minimum of the stationary 
probability distribution, see Fig. 5, e.g.). For this case, according to the 
broad window in the spectrum of eigenvalues found in Sec. 5. (see Fig. 9.), 
the distribution function P(n) relax nearly independently in sections n < ns 
and n > ns after a time defined by second excited eigenvalue ~ l/o~>2, and 
then, after a time of order ~ I/l^i (where uji is first excited eigenvalue) 
the probability start to flow through the unstable point nij, what let the 
distribution function P{n) to relax to the stationary solution Po{n), see 
Fig. 7-8. Motivated by this brief summary of the results presented in Sec. 5., 




-2 2 4 6 8 10 

log t (dimensionless) 



Fig. 10. Average fraction of open channels (n) = / dnnP{n), as a function of time 
for membranes with m = 100 400 channels with symmetric (B = 5, uq = 0.5, 
above) and asymmetric {B = 5, no = 0.49, below) step-like model of ion channel 
switching. The quasi-continuous probability distribution function is defined as 
P(n) = mPk, where n = k/m. The last inflection points on these plots define 
trapping times T'^n), collected in Table 2. 
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we discuss now perturbative approach in terms of weak coupling at the point 
nu- 

Taking into account the probabihty flows between point ku — ^ and 
kjj, where kjj = mriij (see Fig 4), we obtain the equation expressing total 
balance of the probability in the system in the form 

— / P{n, T)dn = nu{l - nu)/^u{T), (29) 

where A[/ is defined as a difference A[/ = P{nu) — P{nu — An), An = 1/m 
and r is an dimensionless time (3). Whereas the cq. (29) has a general 
character, for further discussion the integral Jq^ P{n)dn and difference A[/ 
need to be approximated for r 3> 1/^2 • To establish such an approximation, 
one can define the decomposition of the stationary probability distribution 



-0.002 
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log t (dimensionless) 

Fig. 11. Shanon entropy per channel Ssh/m, where Ssh = / dnP{n)lnP{n), as 
a function of time for membranes with m = 100 400 channels with symmetric 
{B = 5, no = 0.5, above) and asymmetric {B = 5, no = 0.49, below) step-like model 
of ion channel switching. The quasi-continuous probability distribution function 
is defined as P{n) = mPk, where n = k/m. Last inflection points on these plots 
define trapping times Tsh, gathered Table 2. 
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Po = Pl + Pr, where 

p = I Poin) for n<nu p (n) = i ^^^^^ n>nu 

' \ elsewhere ' \ elsewhere 

The representation (30) and a brief discussion of the exact-diagonalization 
results, presented at the beginning of this section, brought us to an ansatz: 

Prjin, t) = Po(n) + riir) [PlH - CflPfl(n)] , (31) 

where ?7(t) is a function of time, that we expect to behave as 77 exp(— wit) 
for r >> l/uJ2, and cr is the coefficient defined to save the normalization 

J^P{n,l)dn = I, CR = Pl/{1 - Pl), where Pl = J^"" Po{n)dn. After 
substituting (31) to the eq. (29), one can write the differential equation for 
the function 77 (r): 

dri nu{l-nu) 

Tr = -P.(l-P.)^^^' 

where Pjj is the value of the stationary probability distribution function 
in point nu, Pjj = Po(n[/). The solution of the eq. (32) is the simple 
exponential decay rj{T) = 77(0) exp(— w^^r), where the decay constant 

_ nu{l-nu) 
"'^ - P.d-PL)"^" ^''^ 

is an approximation of first excited eigenvalue loi. This approximation is 
essentially equivalent to the first-order perturbation correction to the cji, if 
we consider the connection in point no as a perturbation. 

Moreover, one can consider the stationary probability distribution in 
the limit form, for large m (15), and make the Gaussian approximation to 
calculate all the integrals. The practical formula, obtained in this manner 

Ur, = -m{fM - fu) loge + ^ logm + 0{1), (34) 

where /m and fu are values of the function /(n) (see Sec. 4) in point 
um and njj, respectively, can be compared indirectly with numerical exact- 
diagonalization results. Although the coefficient in the linear part of the 
cq. (34) (/m - /t/)loge = 0.01557 and 0.008597, for symmetric (P = 5, 
no = 0.5) and asymmetric {B = 5, no = 0.49) Fermi-Dirac model of the 
ion channel switching, approximately agree with corresponding numerical 
values obtained from exact-diagonalization (see Fig. 9), the term logm/2 
simply does not exists in numerics. This clear contradiction seems to be the 
immanent disadvantage of the perturbation approach to our system, that 
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simply produce too large value of the step A.u in eq. (29). Moreover, it 
cannot be repaired in any perturbation scheme, because, what was shown 
numerically, the second perturbation correction to the first excited eigen- 
value (jJi produce nonphysical negative results. 

To avoid the problem mentioned above, one can derive the diffusion 
equation for n ^ njj: 

dP nu{l - nu) d^P ^^^^ 



dt m dv? ' 



find its Green function 



mm?' 



(36) 



to take the ansatz (31) as the initial condition for the eq. (35) and frozen it 
after a time Tdis- 

Pfi^^ ^diff) = / dn'G{n - n', Tm)Pv{n', r). (37) 

Then, the value of the step Au in eq. (29) is multiplied by factor 2G(0, rjiff) An, 
where G(0, Tdig) is the value of Green's function (36) at point n = and 
An = 1/m. That brought us to the new approximation of first excited 
eigenvalue 



,, froze ^nujl - nu) ApAr, nujl - nu) Pu , . 

= ^^{O^T^dif£)PuAn= d — — — -. (38) 

' Pl{1-Pl) V ^"i^diff Pl{1-Pl) 



The factor appeared in eq. (38) cancel the logarithmic term in (34). 

In so-called frozen diffusion approach, presented above, the diffusion 
time Tdiff is an free parameter, adjusted to agree with constants in for- 
mula for logcji obtained from exact diagonalization. Although the frozen 
diffusion approach produce the results agreeable with the numerics, the 
physical reasons for introducing parameter T^is are not clear at this stage; 
the problem needs a further study, possibly with utilising the Real- Space 
Renormalization Group Technique. 



Summary 

We derived the detail-balance equations for stationary probability dis- 
tribution of the number of open ion channels in the cell membrane. This 
equation was used to compare the extrema of the stationary probability 
distribution for large m-channel system with corresponding stationary so- 
lutions of the deterministic dynamical equation for m — > oo, and to discuss 
the fluctuation around them for finite m. 
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The so-called step-like model of ion channel switching, based on that 
studied by Hodgkin and Huxley [4] , was briefly discussed in relation to the 
experimental data, applied to obtain an analytical solution of the detail- 
balance equation, qualitative discussion of the two classes of phase transi- 
tions in the system, and numerical exact- diagonalization calculations of the 
decay time of the metastable state of the system. The exact diagonalization 
results were shown to disagree with the first-order perturbations. Moreover, 
the perturbation scheme produces unphysical results in the second order. 
This contradiction shows essentially a non-perturbative character of coupled 
ion-channel system. The proposed frozen diffusion approach, with an addi- 
tional adjustable parameter, was shown to explain the exact diagonalization 
results. 
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Appendix A 

Fluctuations around the stationary points obtained from the 
detailed-balance approach 

In this appendix, we present the approximate derivation of the small 
fluctuations of the number of open channels around the stationary points 
for general model of ion channel switching. Let us apply the detailed balance 
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equation (13) to calculate the ratio of the probability distribution Pk^/^/Pk 
for points in distance A: 

D A-l J-, A-i , • 1 

Pk-A _ JT Pk-j-l _ JT k-J 1 - Pk-j 

Pk f}^ Pk-j /j, m-k + j + l pk-,-1 ' ^ ■ ^ 

where probability of open channel pk = Popenik/m). For A/m « 1 (small 
fluctuations), which supposition would be verified a posteriori, one can ap- 
proximate the probabilities in eq. (A.l) 

Pk-j ^Pk- jp'k = Pfe(l - aj), 

where = p'openij 1"^) /f^^ and 

a = ^ = (A.2) 
Pk rnpopen{k/m) 

and analogically 

1 - Pk-j ~ 1 - Pjk(l - aj) = (1 - Pjfc)(l + Pj), 

where: 

/? = a-^ = -^«l. (A.3) 

l-Pk l-Pk 

Subsequently, eq. (A.l) can be rewrite in the approximate form: 
Pk-A^, kl (m-kV. fl-Pk)\\^^^ 



Pk {k-A)\im-k + A)\\ Pk 
where 11 a is the non-binomial factor in ratio Pk-A/Pk- 



"A = n 1 _ Z-i u " n [1 + - - 



(a+/3)AV2 



(A.5) 



where we skipped the terms of order A/m. 

Furthermore, for stationary point pk = k/m = n^, (see Sec. 3). So that, 
using Stirling approximation In/c! ^ kink — k, one can obtain from eq. 
(A.4): 

Pk-A 

In—— ^ -MuHa. 

Pk 2mn*(l — n*) 

After substituting eq. (A.5) and, subsequently (A.2) and (A.3), we finally 
obtain: 

In « (A.6) 

Pk 2mn*(l — n*) 
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where p'^ = p'openi''^*)- -^l- (-^-6) let us read the mean-square fluctuation of 
the fraction of open channels 6n = A/m around the stationary point n*: 

\ / m 1 — p'^ 

that confirms our starting supposition A/m << 1 for large m. Moreover, eq. 
(A.7) let us to identify the value p'^ = 1, previously distinguishing between 
types of stationary sohitions (see Sec. 3.), with critical point in terms of 
2-nd order phase transitions, for which 5n ^ oo. Additionally, for quasi- 
continuous probability distribution function P{n), defined by the relation 
Pfe = P{n)An (cf. Sec. 4.) one can write (for n n*): 



P(n) = P(n*) exp 



—m 



1-K 
n*(l — 



-(n - n*)^, 



(A.8) 



what clearly identifies stable and unstable solutions of the relation = 
Popen{iT'*) with, respectively maxima and minima of the probability distri- 
bution. Eq. (A.8) motivates the ansatz 

P ^ e^fin) in Sec. 4. 



